# Markdown
library(bookdown)
# General data analysis and transformation
library(readr)
library(readxl)
library(dplyr)
library(janitor)
# library(gtools)
library(stringr)
library(units)
library(tidyr)
library(rlang)
# API-related
library(jsonlite)
library(urltools)
# Mapping and GIS operations
library(sf)
library(leaflet)
library(htmlwidgets)
library(geojsonsf)
library(rmapshaper)
library(tmap)
library(tmaptools)
library(ceramic)
# workflow
library(here)
# plotting
library(ggplot2)
library(forcats)
# Get list of CES 3 variable names and plain text names
ces3_variables <- readr::read_csv(here('data_prepared', 'ces_names.csv'))
ces_measures <- ces3_variables %>%
filter(group == 'pollution burden',
subgroup == 'exposures',
type == 'percentile')
This document attempts to do some quantitative geospatial analysis to look at correlations between redline mapping and indicators of water quality, public health, and facility information. It computes the area weighed average percentile for different CalEnviroScreen 3.0 (CES) indicators by HOLC rating for a given city and statewide. This document considers indicators in the ‘Exposures’ category of the CES Pollution Burden indicators (factsheet, indicators page), including:
The CES polygons and the Redline polygons have different coverages, as shown in the maps below. In the first map for each indicator, the left pane of each map shows the HOLC rated areas, the center pane shows the CES indicator scores (at the census tract level) with the Redlined areas superimposed on top, and the right pane shows the overlapping portions of the CES polygons and the redlined areas (border colors represent HOLC rating).
# load data
# Redline Polygons
redline_polygons <- read_rds(here('data_prepared', 'redline_polygons.RDS'))
# st_crs(redline_polygons) # to check the reference system
# modify the 'city' column to better keep track when joining with other datasets
redline_polygons <- redline_polygons %>% rename('redline_city' = 'city')
# add the area of each redline polygon
redline_polygons <- redline_polygons %>% mutate(redline_poly_area = st_area(.))
# CES 3 Polygons
ces_3_poly <- read_sf(here('data_raw', 'CalEnviroScreen3', 'CESJune2018Update_SHP','CES3June2018Update.shp'))
# st_crs(ces_3_poly) # to check the reference system
# modify the 'City' column to better keep track when joining with other datasets
ces_3_poly <- ces_3_poly %>% rename('CES_City' = 'City')
# Transform both to UTM Zone 10N
redline_polygons <- redline_polygons %>% st_transform(crs = 26910)
ces_3_poly <- ces_3_poly %>% st_transform(crs = 26910)
# create a map to check that the data loaded correctly
# tm_shape(ces_3_poly) + tm_borders() + tm_shape(redline_polygons) + tm_borders(col = 'red')
# get a basemap (See: https://github.com/mtennekes/tmap/issues/185)
Sys.setenv(MAPBOX_API_KEY = "pk.eyJ1IjoiZGFsdGFyZSIsImEiOiJjazZqcm5ocHgwMG84M2tteGc5eWd4YXNmIn0.nHiYenoq-wFWwH3V4vtPXA")
background <- cc_location(ces_3_poly %>% filter(CES_City == 'Sacramento'), verbose = FALSE)
# create map
map_redline <- tm_shape(background) +
tm_rgb(alpha = 0.5) +
tm_shape(redline_polygons %>% filter(redline_city == 'Sacramento'), is.master = TRUE) +
tm_polygons('holc_grade', palette = c('green', 'blue', 'yellow', 'red'), title = 'HOLC Grade') +
tm_layout(legend.bg.color = 'white', legend.bg.alpha = 0.7, legend.frame = TRUE)
# get the intersection of the CES polygons and redline polygons
ces3_redline_intersection <- st_intersection(ces_3_poly, redline_polygons)
# inspect
# glimpse(ces3_redline_intersection)
# add the area of each clipped polygon to the data frame
ces3_redline_intersection <- ces3_redline_intersection %>%
mutate(clipped_area = st_area(.))
# glimpse(ces3_redline_intersection)
ces_poly_map <- function(ces_id, ces_measure_name) {
# plot the CES polygons
map_ces <- tm_shape(background) +
tm_rgb(alpha = 0.5) +
tm_shape(ces_3_poly %>% filter(CES_City == 'Sacramento' | CES_City == 'West Sacramento')) +
tm_polygons(col = ces_id, title = ces_measure_name, border.alpha = 1) +
tm_shape(ces3_redline_intersection %>% filter(redline_city == 'Sacramento'), is.master = TRUE) +
tm_borders(lwd = 0, alpha = 0) +
# tm_shape(redline_polygons %>% filter(redline_city == 'Sacramento')) +
# tm_borders(lwd = 2, col = 'black') +
tm_shape(redline_polygons %>%
filter(redline_city == 'Sacramento', holc_grade == 'D')) +
tm_borders(lwd = 2, col = 'red') +
tm_shape(redline_polygons %>%
filter(redline_city == 'Sacramento', holc_grade == 'C')) +
tm_borders(lwd = 2, col = 'gold2') +
tm_shape(redline_polygons %>%
filter(redline_city == 'Sacramento', holc_grade == 'B')) +
tm_borders(lwd = 2, col = 'blue') +
tm_shape(redline_polygons %>%
filter(redline_city == 'Sacramento', holc_grade == 'A')) +
tm_borders(lwd = 2, col = 'green') +
# tm_text('holc_grade', size = 0.5) +
tm_layout(legend.bg.color = 'white', legend.bg.alpha = 0.7, legend.frame = TRUE)
map_ces
}
ces_redline_overlap_map <- function(ces_id, ces_measure_name) {
map_overlap <- tm_shape(background) +
tm_rgb(alpha = 0.5) +
tm_shape(ces_3_poly %>% filter(CES_City == 'Sacramento')) +
tm_borders(lwd = 1) +
tm_shape(ces3_redline_intersection %>% filter(redline_city == 'Sacramento'), is.master = TRUE) +
tm_polygons(col = ces_id, title = ces_measure_name, border.alpha = 0) +
# tm_shape(redline_polygons %>% filter(redline_city == 'Sacramento')) +
# tm_borders(lwd = 2, col = 'black') +
tm_shape(redline_polygons %>%
filter(redline_city == 'Sacramento', holc_grade == 'D')) +
tm_borders(lwd = 2, col = 'red') +
tm_shape(redline_polygons %>%
filter(redline_city == 'Sacramento', holc_grade == 'C')) +
tm_borders(lwd = 2, col = 'gold2') +
tm_shape(redline_polygons %>%
filter(redline_city == 'Sacramento', holc_grade == 'B')) +
tm_borders(lwd = 2, col = 'blue') +
tm_shape(redline_polygons %>%
filter(redline_city == 'Sacramento', holc_grade == 'A')) +
tm_borders(lwd = 2, col = 'green') +
#tm_text('holc_grade', size = 0.5) +
tm_layout(legend.bg.color = 'white', legend.bg.alpha = 0.7, legend.frame = TRUE)
map_overlap
}
# FACET MAP
ces_redline_facet_map <- function(ces_id, ces_measure_name) {
map_facets <- tm_shape(ces_3_poly %>% filter(CES_City == 'Sacramento')) +
tm_borders(lwd = 1) +
tm_shape(background) +
tm_rgb(alpha = 0.5) +
tm_shape(ces3_redline_intersection %>%
filter(redline_city == 'Sacramento'), is.master = TRUE) +
# tm_shape(ces3_redline_intersection %>% filter(redline_city == 'Sacramento')) +
tm_facets(by = 'holc_grade', free.coords = FALSE) +
tm_polygons(col = ces_id, title = ces_measure_name, border.alpha = 0) +
tm_shape(redline_polygons %>% filter(redline_city == 'Sacramento')) +
tm_facets(by = 'holc_grade') +
tm_borders(lwd = 2, col = 'black') +
# # D
# tm_shape(redline_polygons %>%
# filter(redline_city == 'Sacramento', holc_grade == 'D')) +
# tm_facets(by = 'holc_grade') +
# tm_borders(lwd = 2, col = 'red') +
# # C
# tm_shape(redline_polygons %>%
# filter(redline_city == 'Sacramento', holc_grade == 'C')) +
# tm_borders(lwd = 2, col = 'gold2') +
# tm_facets(by = 'holc_grade') +
# # B
# tm_shape(redline_polygons %>%
# filter(redline_city == 'Sacramento', holc_grade == 'B')) +
# tm_borders(lwd = 2, col = 'blue') +
# tm_facets(by = 'holc_grade') +
# # A
# tm_shape(redline_polygons %>%
# filter(redline_city == 'Sacramento', holc_grade == 'A')) +
# tm_borders(lwd = 2, col = 'green') +
# tm_facets(by = 'holc_grade') +
#
tm_legend(bg.color= 'white', frame = TRUE)
map_facets
}
ces_scores_plot <- function(ces_id, ces_measure_name) {
# group by city and holc rating, then compute the weighted average score for each
ces3_redline_grouped <- ces3_redline_intersection %>%
st_drop_geometry() %>%
mutate(area_x_score = clipped_area * (!!as.name(ces_id))) %>%
group_by(redline_city, holc_grade) %>%
summarize(total_area = sum(clipped_area),
area_x_score_total = sum(area_x_score)) %>%
mutate(weighted_score = area_x_score_total / total_area) %>%
drop_units()
ces3_redline_grouped_statewide <- ces3_redline_grouped %>%
group_by(holc_grade) %>%
summarize(total_area = sum(total_area),
area_x_score_total = sum(area_x_score_total)) %>%
mutate(weighted_score = area_x_score_total / total_area) %>%
mutate(redline_city = '**Statewide**') %>%
drop_units()
ces3_redline_grouped_city_state <- bind_rows(ces3_redline_grouped, ces3_redline_grouped_statewide)
# plot the results by city and holc rating, and include statewide summary
g_city_state <- ggplot(ces3_redline_grouped_city_state) +
aes(x = fct_relevel(fct_reorder(redline_city, weighted_score),
c('**Statewide**'),
after = 0),
y = weighted_score,
fill = fct_rev(holc_grade)) +
geom_bar(stat = 'identity', position = 'dodge') +
geom_vline(xintercept = 1.52, size = 0.5, color = 'grey50') +
scale_fill_manual(values = rev(c('green', 'blue', 'yellow', 'red'))) +
guides(fill = guide_legend(reverse = TRUE)) +
labs(fill = 'HOLC Rating',
title = ces_measure_name,
x = NULL,
y = "Area Weighted Average Score By City") +
coord_flip()
g_city_state
}
i_map <- 1
ces_measure_name <- ces_measures$name[i_map]
ces_id <- ces_measures$id[i_map]
map_ces <- ces_poly_map(ces_id, ces_measure_name)
map_overlap <- ces_redline_overlap_map(ces_id, ces_measure_name)
tmap_arrange(map_redline, map_ces, map_overlap, nrow = 1)
ces_redline_facet_map(ces_id, ces_measure_name)
ces_scores_plot(ces_id, ces_measure_name)
i_map <- 2
ces_measure_name <- ces_measures$name[i_map]
ces_id <- ces_measures$id[i_map]
map_ces <- ces_poly_map(ces_id, ces_measure_name)
map_overlap <- ces_redline_overlap_map(ces_id, ces_measure_name)
tmap_arrange(map_redline, map_ces, map_overlap, nrow = 1)
ces_redline_facet_map(ces_id, ces_measure_name)
ces_scores_plot(ces_id, ces_measure_name)
i_map <- 3
ces_measure_name <- ces_measures$name[i_map]
ces_id <- ces_measures$id[i_map]
map_ces <- ces_poly_map(ces_id, ces_measure_name)
map_overlap <- ces_redline_overlap_map(ces_id, ces_measure_name)
tmap_arrange(map_redline, map_ces, map_overlap, nrow = 1)
ces_redline_facet_map(ces_id, ces_measure_name)
ces_scores_plot(ces_id, ces_measure_name)
i_map <- 4
ces_measure_name <- ces_measures$name[i_map]
ces_id <- ces_measures$id[i_map]
map_ces <- ces_poly_map(ces_id, ces_measure_name)
map_overlap <- ces_redline_overlap_map(ces_id, ces_measure_name)
tmap_arrange(map_redline, map_ces, map_overlap, nrow = 1)
ces_redline_facet_map(ces_id, ces_measure_name)
ces_scores_plot(ces_id, ces_measure_name)
i_map <- 5
ces_measure_name <- ces_measures$name[i_map]
ces_id <- ces_measures$id[i_map]
map_ces <- ces_poly_map(ces_id, ces_measure_name)
map_overlap <- ces_redline_overlap_map(ces_id, ces_measure_name)
tmap_arrange(map_redline, map_ces, map_overlap, nrow = 1)
ces_redline_facet_map(ces_id, ces_measure_name)
ces_scores_plot(ces_id, ces_measure_name)
i_map <- 6
ces_measure_name <- ces_measures$name[i_map]
ces_id <- ces_measures$id[i_map]
map_ces <- ces_poly_map(ces_id, ces_measure_name)
map_overlap <- ces_redline_overlap_map(ces_id, ces_measure_name)
tmap_arrange(map_redline, map_ces, map_overlap, nrow = 1)
ces_redline_facet_map(ces_id, ces_measure_name)
ces_scores_plot(ces_id, ces_measure_name)
i_map <- 7
ces_measure_name <- ces_measures$name[i_map]
ces_id <- ces_measures$id[i_map]
map_ces <- ces_poly_map(ces_id, ces_measure_name)
map_overlap <- ces_redline_overlap_map(ces_id, ces_measure_name)
tmap_arrange(map_redline, map_ces, map_overlap, nrow = 1)
ces_redline_facet_map(ces_id, ces_measure_name)
ces_scores_plot(ces_id, ces_measure_name)